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The liquid-liquid critical point scenario of water hypothesizes the existence of two metastable liq- 
uid phases — low-density liquid (LDL) and high-density liquid (HDL) — deep within the supercooled 
region. The hypothesis originates from computer simulations of the ST2 water model, but the sta- 
bility of the LDL phase with respect to the crystal is still being debated. We simulate supercooled 
ST2 water at constant pressure, constant temperature and constant number of molecules A'' for 
N < 729 and times up to 1 fis. We observe clear differences between the two liquids, both structural 
and dynamical. Using several methods, including finite-size scaling, we confirm the presence of a 
liquid-liquid phase transition ending in a critical point. We find that the LDL is stable with respect 
to the crystal in 98% of our runs (we perform 372 runs for LDL or LDL-like states), and in 100% 
of our runs for the two largest system sizes (A^ = 512 and 729, for which we perform 136 runs 
for LDL or LDL-like states). In all these runs tiny crystallites grow and then melt within 1 fis. 
Only for N < 343 we observe six events (over 236 runs for LDL or LDL-like states) of spontaneous 
crystallization after crystallites reach an estimated critical size of about 70 ± 10 molecules. 

PACS numbers: 64.60. F-, 64.70.Ja, 82.60.-s, 07.05.Tp, 61.20.Ja, 61.25.Em 



I. INTRODUCTION 

For many centuries, water and its anomalies have been 
of much interest to scientists. A particular rise of inter- 
est occurred in the late 1970s after experiments done by 
Angell and Speedy seemed to imply some kind of criti- 
cal phenomenon in supercooled liquid water at very low 
temperatures [IHlj. Even though liquid water experi- 
ments are limited by spontaneous crystallization below 
the homogenous nucleation temperature {Th ~ 233 K 
at 1 bar), it is possible to further explore the phase dia- 
gram by quenching water to far lower temperatures [5HI] ■ 
The result of these experiments is an amorphous solid, 
i.e. a glassy ice, corresponding to an out-of equilibrium 
state that is very stable with respect to the equilibrium 
crystalline ice phase. The amorphous depends on the ap- 
plied pressure: at low pressure, below w 0.2 GPa, the low 
density amorphous ice (LDA) is formed, while at higher 
pressure the high density amorphous ice (HDA) is ob- 
served [S]. It has been shown by Mishima et at that 
these two amorphous ices are separated by a reversible 
abrupt change in density that resembles in all its respects 
an equilibrium first order phase transition [941 2j. 

Raising the temperature of either LDA or HDA does 
not turn the sample into a liquid, but leads once again 
to spontaneous crystallization (around Tx ~ 150 K). In 
fact, between and Tx, often called the "no man's 
land" of bulk water, crystallization occurs at a time 
scale that is too short for current experimental meth- 
ods, although a new technique is possibly succeeding in 
the task of measuring the metastable liquid phase |13j . 



Computer simulations of water, however, involve time 
scales small enough to witness spontaneous crystalliza- 
tion and are therefore able to explore liquid water in the 
"no man's land". In 1992 Poole et al. [HI performed a 
series of molecular dynamics simulations using the ST2 
water model [15], using the reaction field method for 
the long-range interactions (ST2-RF), and discovered a 
liquid-liquid phase transition ending in a critical point, 
separating a low density liquid (LDL) and a high den- 
sity liquid (HDL). These two liquids can be considered 
to be the liquid counterparts of the LDA and HDA, re- 
spectively. 

The existence of the critical point also allows one to 
understand X-ray spectroscopy results [rMTO] . explains 
the increasing correlation length in bulk water upon cool- 
ing as found experimentally ^20], the hysteresis effects 
[21) and the dynamic behavior of protein hydration water 
[^^24| . It would be consistent with a range of thermo- 
dynamical and dynamical anomalies [25H33j and experi- 
ments j34-37 . 

Many more computer simulations investigating the 
phenomenology of the liquid-liquid critical point (LLCP) 
have been performed since then j55H55] . Detailed studies 
using ST2-RF have been made by Poole et al. [SB] using 
molecular dynamics, while Liu et al. simulated ST2 with 
Ewald summation (ST2-Ew) for the electrostatic long- 
range potential using Monte Carlo [57] [35] • Also in other 
water models the liquid-liquid phase transition (LLPT) 
and its LLCP are believed to be found, for example by 
Yamada et al. in the TIP5P model [59], by Paschek et 
al. in the TIP4P-Ew model [SO], and in TIP4P/2005 by 



Abascal and Vega [5T1 . 

Recently Limmer and Chandler used Monte Carlo um- 
brella sampling to investigate the ST2-Ew model, but 
claimed to have found only one liquid metastable phase 
(HDL) rather than two [63j. They therefore concluded 
that LDL does not exist because it is unstable with re- 
spect to either the crystal or the HDL phase. The em- 
phasis in their work is about the difference between a 
metastable phase, i.e. separated from the stable phase by 
a finite free-energy barrier, and an unstable state, where 
the free-energy barrier is absent and the state does not 
belong to a different phase. 

Shortly after, Poole et al. [M] and Kesselring et al. 
[65j presented results using standard molecular dynam- 
ics for ST2-RF showing the occurrence of the LLCP with 
both HDL and LDL phases metastable with respect to 
the crystal, but with the LDL not unstable with respect 
to either the crystal or the HDL. This result was con- 
firmed, using the same method as Limmer and Chandler, 
by Sciortino et al. [66 and Poole et al. 67J in ST2-RF 
water and by Liu el al. in ST2-Ew water [68J. 

The aim of this paper is to confirm the presence of a 
liquid-liquid critical point in water in the thermodynamic 
limit using finite size scaling techniques, and confirm that 
LDL is a bona fide metastable liquid. We use the ST2- 
RF model because it has been well-studied in the super- 
cooled region, making it easier to compare and verify our 
data. In the supercooled phase it has a relatively large 
self-diffusion compared to other water models, therefore 
suffers less from the slowing down of the dynamics at ex- 
tremely low temperatures. We explore a large region of 
the phase diagram of supercooled liquid ST2-RF water 
(Fig. [l]) using molecular dynamics simulations with four 
different system sizes by keeping constant the number 
N of molecules, the pressure P and the temperature T 
{NPT ensemble). 

Within the explored region we find both LDL and 
HDL, separated at high pressures by a LLPT, ending 
in a LLCP estimated at Pc w 208 MPa and Tc ~ 246 K. 
This phase transition is particularly clear in Fig. [2] where 
one can see from the density how the system continu- 
ously flips between the two states. However, due to finite 
size effects this phase flipping also occurs below the crit- 
ical point along the Widom line (the locus of correlation 
length maxima) [SHI [7D] . For this reason it is necessary 
to apply finite size scaling methods to establish the exact 
location of the critical point. 

For six state points and for small system size N < 343 
we observe, in only one over the (on average) seven sim- 
ulations we performed for each state point, irreversible 
crystal growth, indicated as red stars in Fig. [T] Each 
of these crystallization events occurred within the LDL 
(or LDL-like) region. Analysis of these crystals revealed 
them to have a diamond cubic crystal structure. As we 
will discuss later, because these events disappears for 
larger systems, we ascribe these crystallization to finite- 
size effects. 

We start in Sec. IH] with a description of the model 
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FIG. 1: Overview of the state points at which simulations 
have been performed, with the symbols indicating different 
system sizes. At high temperatures we observe a high-density 
liquid state (HDL, shaded in pink), while at lower temper- 
atures we find a low-density liquid (LDL, in blue). These 
are separated by a region where the system is continuously 
flipping between the two states, as seen in Fig. [2] This tran- 
sition region (in purple) is identified as the liquid-liquid phase 
transition line (LLPT) at high pressures, and the Widom 
line at low pressures. These lines join at the liquid-liquid 
critical point (LLCP) estimated at Pc = 208 ± 3 MPa and 
Tc = 246± 1 K (see Sec. |VII[ ). At low temperatures the LDL 
(or LDL-like) region is bound by the glass transition temper- 
ature Tg, below which we can no longer fully equilibrate the 
system within 100 ns, and consider the liquid to have become 
a glass (see Sec. IV I. For small sizes (A'^ < 343) we observe 
spontaneous crystallization within 1 ns-long simulations at six 
state points (indicated by the red stars), all of them within 
the LDL (or LDL-like) region. We never observe crystalliza- 
tion for sizes A'^ = 512 and 729 for simulations of comparable 
duration. Because the probability of crystalization should 
increase with N, this results suggest that our cystallization 
events are a finite-size effect that becomes negligeble for large 
sizes. Crystallization events are discussed in Sec|VI| 



HI we dis- 



and the procedures that were used. In Sec. 
cuss the use of the intermediate scattering function to 
analyze the structure of the liquid, and in Sec. |IV| its 
use in defining the correlation time. The analysis of the 
liquid structure is continued in Sec. |V] where we define 
and compare a selection of structural parameters. The 
parameter is found to be particularly well-suited to 
distinguish between the liquid and the crystal state, and 
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FIG. 2: Phase flipping near the phase transition line (P = 
215 MPa, with N = 343 molecules). At high T the sys- 
tem is in the HDL phase (with a density p ~ 1.03 g/cm''), 
while at low T the system is in the LDL phase (density 
p ~ 0.88 g/cm^). However, near the phase transition line (at 
T ~ 244.5 K for this pressure) the system is flipping between 
the two phases. 



this fact is subsequently used in Sec. |VI| where we dis- 
cuss the growth and melting of crystals within the LDL 
liquid. In Sec. |VII[ by defining the appropriate order pa- 
rameter, we show that the LLCP in ST2-RF belongs to 
the same universality class as the 3D Ising model. We 
accurately determine where the LLCP is located in the 
phase diagram in the thermodynamic limit by applying 
finite size scaling on the Challa-Landau-Binder parame- 
ter. We discuss our results and present our conclusions 
in Sec. Ivml 



II. SIMULATION DETAILS 

In the ST2 model [H] each water molecule is rep- 
resented by a rigid tetrahedral structure of five parti- 
cles. The central particle carries no charge and rep- 
resents the oxygen atom of water. It interacts with 
all other oxygen atoms via a Lennard-Jones (LJ) po- 
tential, Ui,i{rij) EE 4£ [(cr/r,j)i2 - {a/r,jf] with e = 
0.31694 kJ/mol and a = 3.10 A. Two of the outer parti- 
cles represent the hydrogen atoms. Each of them carries 
a charge of -1-0.2357 e, and is located a distance 1 A away 
from the central oxygen atom. The two remaining parti- 
cles carry a negative charge of —0.2357 e, are positioned 
0.8 A from the oxygen, and represent the lone pairs of a 
water molecule. 

The electrostatic potential in ST2 is treated in a special 
way. To prevent charges a and b from overlapping, the 
Coulomb potential is reduced to zero at small distances: 



C^cl(''afc) = S{rij) 
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where S{rij) is a function that smoothly changes from 
one to zero as the distance between the molecules de- 
creases. 
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with Rl = 2.0160 A, Ru = 3.1287 A, and where r,j is 
the distance between the oxygen atoms of the interacting 
molecules. In the original model a simple cutoff was used 
for the electrostatic interactions. In this paper, however, 
we apply the reaction field method |71) which changes 
the ST2 Coulomb potential to 

^el(...)^5(..,m..)|^(^ + ^) (3) 

where T{rij) is another smoothing function: 
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We use a reaction field cutoff Rc = 7.8 A together with 
Rt = 0.95i?c. These parameters define our ST2-RF wa- 
ter model and are the same that were used in previous 
ST2-RF simulations. 

For the LJ interaction we use a simple cutoff at the 
same distance of 7.8 A. We do not adjust the pressure to 
correct for the effects of the L J cutoff [72l [73] , since these 
adjustments come from mean field calculations which be- 
come increasingly weak as one approaches a critical point. 

We use the SHAKE algorithm [74| to keep the relative 
position of each particle within a ST2 molecule fixed. The 
temperature and pressure are held constant using a Nose- 
Hoover thermostat [731 [731 [75] together with a Berendsen 
barostat [77] • In all simulations periodic boundary con- 
ditions are applied. 

Our code is validated by simulating the same state 
points as those published by Poole et ai, see Fig. lb in 
[56], where pressure corrections for the LJ cutoff were 
applied in the NVT (constant iV, T and volume V) en- 
semble. Averaging at each state point over 10 simulations 
with different initial conditions allows us to estimate the 
error bars. In Fig. [Sjwe compare our results for N — 216 
molecules and density 0.83 g/av? , and find that our data, 
after pressure correction, matches that of Ref. [56 well. 

For each of the simulations done in the NPT ensem- 
ble, we use the following protocol. We first create a box 
of N molecules at n different initial densities (with n up 
to 21) ranging from 0.85 to 1.05 g/cm^ . We then perform 
a 1 ns NVT simulation at T = 300 K. In this way we 
obtain n independent configurations all &t T — 300 K 
in the prefixed range of densities. Next, we use these 
independent configurations as starting points for NPT 



4 



-70 r 



-80- 



-90- 



1 -100- 



-110- 




■ Poole, el al. (Fig. lb) 

O NVT incl. pressure corrections 



-120 L 



_1_ 



260 280 300 320 
Temperature T [K] 



340 



360 



FIG. 3: To validate our code, we compare our simula- 
tion results with those from Poole et al. [SB] at density p = 
0.83 g/cva? and for A'^ = 216 molecules. We performed simual- 
tions in the NVT ensemble applying pressure corrections and 
find the same results as Ref. [56] within the error bars. At this 
density the pressure correction due to the LJ cutoff (propor- 
tional to p^) is equal to —12.66 MPa. The variation of P with 
T along this isochore shows the occurence of both a density 
maximum at 300 K and a density minimum near 265 K, as 
at these state points {dp/dT)p = -pKT{dP/dT)v = with 
Kt > the isothermal compressibility. 



simulations at T = 265 K and different pressures ranging 
from 190 to 240 MPa, and continue the simulation for 
an additional 1 ns. This results in n independent con- 
figurations at T = 265 K and the given pressure. For 
all pressures considered here, this will lead the system 
into the HDL phase. Finally the system is quenched to 
the desired temperature at the given pressure, followed 
by 100-200 ns of equilibration time. 



In Sec. |IV| it will 
be shown that this provides enough time for the system 
to reach equilibrium for the state points above the line 
marked with the label Tg in Fig. [T] 



III. INTERMEDIATE SCATTERING 
FUNCTION 



The intermediate scattering function ^(k, t) plays an 
essential role in the analysis of liquid structure, since it 
is frequently measured in experiments as well as easily 
calculated from simulation data. It describes the time 
evolution of the spatial correlation at the wave vector k, 
and can be used to distinguish between phases of different 
structure, such as LDL and HDL or crystal. It is defined 
as 

5(k,i) = - /^e'k.[r,(t')-r™(*'+t)]\ 



where denotes averaging over simulation time t', 

and vi{t') the position of particle t at time t' . For sim- 



plicity we only apply the intermediate scattering function 
to the oxygen atoms, which we denote as Soo{k,i). 

Since the system has periodic boundary conditions, the 
components of k have discrete values 27rn/L, where L 
is the length of the simulation box and n = 1,2,3,.... 
We define Sooik,t) = (S'oolk, i))„ where the average 
is taken over all vectors k with magnitude k belonging 
to the nth spherical bin n{n — \)/ L < k < TT{n + \) / L 
for n — 2,3, .. . ,300. Similarly, we define the structure 
factor Soo{k) = {Sooik,t))^ as the time-averaged inter- 
mediate scattering function, with (unless indicated oth- 
erwise) the average taken over the whole duration of the 
run. 
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FIG. 4: The structure factor Soo{k) for a range of temper- 
atures at (a) 210 MPa and (c) 200 MPa for iV = 729. (a) For 
P > Pc the structure has a large change between T = 245 
and 246 K, corresponding to the LDL-HDL first-order phase 
transition, (b) The value of Soo for k corresponding to the 
first maximum, the first minimum and the second maximum 
as a function of T for P = 210 MPa as in panel (a), (c) For 
P < Pc the structure changes in a way that is smoother than 
the case in panel (a), with the more evident change occurring 
between T = 249 and 250 K, corresponding to the crossing 
of the Widom line, as marked by the value of Soo at first 
maxima and minima in panel (d). 

We study Soo{k) above and below our estimate for the 
LLCP pressure. At P = 210 MPa > Pc (Fig. |4|i,b) we 
observe a discontinuous change in the first two peaks of 
Soo{k) as T changes between 245 and 246 K, and a con- 
tinuous change above and below this temperatures. This 
is the expected behavior for a first order phase transi- 
tion occurring at 245 K< T < 246 K and P = 210 MPa 
between two phases with different structure, consistent 
with our results in Fig. [l] The fact that for both phases 
Soo{k) ^ 0{1) for all k shows that both phases are fiuid. 
Indeed, for a crystal-like configuration, with a long-range 
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order, there would be at least one wave vector such that 
Soo{k) ~ 0{N) [78]. Furthermore, the fact that at lower 
T the first peak increases and the other peaks only have 
minor changes indicates that the lower-T liquid has a 
smaller density than the highcr-T liquid. Therefore, this 
result show a first-order phase transition between the 
LDL at lower-T and HDL at higher-T. This transition 
occurs at the same temperature at which we observe the 
phase flipping in density (Fig. [2]) and corresponds to the 
purple region at P > Pc in Fig. [TJ 

The fact that the peaks of Soo{k) are sharper in LDL 
than HDL is an indication that the LDL phase is more 
structured. We can also observe that the major struc- 
tural changes in Sooik) between LDL and HDL are for 
k ~ 1.8 and 2.8 A~^, corresponding to r = in/k ~ 7 and 
4.5 A, respectively, i.e. are for the third and the second 
neighbor water molecules. This change in the structure is 
consistent with a marked shift inwards of the second shell 
of water with increased density, and almost no change in 
the first shell (at k ~ 4.6 and r ~ 2.75 A), as seen in 
structural experimental data for supercooled heavy wa- 
ter interpreted with Reverse Monte Carlo method [75] . 
This changes are visible also in the 00 radial distribu- 
tion function gooir) (Fig. [5^,b). 
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transforming continuously in a shoulder. Same qualita- 
tive behavior is observed for gooi^) (Fig- [sj^id). These 
quantities show us also that the lower-T structure is LDL- 
like, while the higher-T structure is HDL-like. However, 
the absence of any discontinuous change in the struc- 
ture implies the absence of a first-order phase transition 
in the structure of the liquid. This is consistent with 
the occurrence of a LLCP at the end of the first-order 
phase transition somewhere between 200 and 210 MPa, 
at a temperature between 245 and 250 K. In Sec. |VII| we 
shall apply a different method to locate the LLCP with 
more precision. 

At P < Pc, m the one-phase region, we expect to find 
the Widom line emanating from the LLCP. The Widom 
line is by definition the locus of maxima of the correla- 
tion length, therefore, for general thermodynamic con- 
siderations [70] near the LLCP it must be also the lo- 
cus of maxima of the response functions. In particular, 
it must be the locus where the isobaric heat capacity 
Cp = T{dS/dT)p, where S is the entropy of the sys- 
tem, has its maximum along a constant-P path. This 
maximum occurs where the entropy variation with T is 
maximum, expected where the structural variation of the 
liquid is maximum, i.e. where the derivatives of the val- 
ues of Soo{k) (Fig.|4|l) and goo{r) (Fig.[5]i) with T are 
maximum. The interval of temperatures for each P where 
this occurs corresponds to the purple region at P < Pc 
in Fig. [T] indicated as the Widom line. 

It is actually possible to follow the structural changes 
during the simulation. An example is given in Fig. [6] 
where we focus on a 30 ns time period of a simulation at 
200 MPa and 248 K. We divide this time period into six 
5 ns intervals and for each interval we calculate the in- 
termediate scattering function, time-averaged over those 
5 ns. We observe that the liquid is LDL-like for the 
first and third interval, having low density and LDL-like 
Soo{k) (first peak near 2 A~^, separated from the sec- 
ond). On the contrary, for the fifth and sixth interval the 
density is high and Soo{k) is HDL-like (the first peak is 
merely a shoulder of the second peak), indicating that 
the liquid is HDL-like. For the second and fourth in- 
terval, the liquid has an intermediate values of density 
and Sooik), indicating that it is a mix of LDL-like and 
HDL-like structures. 



FIG. 5: The radial distribution function gooif) for the state 
points in Fig. |4] (a) For P > Pc the main structural change 
between LDL and HDL is visible around the second coordina- 
tion shell at r ~ 4.6 A and is stronger when T changes between 
245 and 246 K, as emphasized by the change of values of goo 
for r corresponding to the first maximum and minimum and 
the second maximum in panel (b). (c) The transition from 
LDL to HDL is smoother for P < Pc when the system is 
crossing the Widom line, as shown by the variation of the 
values of goo in panel (d). 

For P < Pc (Fig. |4]:,d) by increasing T we observe 
that the first peak of Soo{k) merges with the second. 



IV. CORRELATION TIME 

Apart from its use in structure analysis, the interme- 
diate scattering function Soo{k,t) can also be used to 
define a correlation time r, i.e. the time it takes for a 
system to lose most of its memory about its initial con- 
figuration 80, 81J. 

In Fig. [t] we show how Sooik, t) decays with time 
for a fixed value of k. Its decay is characterized by 
two relaxation times, the a-relaxation time and the 
/3-relaxation time r^. On very short time scales, the 
molecules do not move around much and each molecule is 
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FIG. 6: As the density changes from p(LDL) to p(HDL), 
also the structure changes. The inset shows how the density 
is changing with time for six consecutive time intervals of 
10 ns, with the corresponding Soo{k) shown in the main plot 
(N = 343 at 200 MPa and 248 K). 
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FIG. 7: Decay of Soo{k,t) with time, for P = 210 MPa, 
T = 250 K and N = 343. Symbols indicate Foo{ki,t) for 
three different values of k: the first maximum of Sooik) at 
ki (red circles), the second maximum at k2 (blue squares), 
and the third maximum fca (green diamonds). Solid lines are 
fits according to Eq. ([sjl. The two components of Eq. ([5| are 
explicitly shown for Foo{k3,t): the green dashed line repre- 
sents the /3-relaxation and is given by [1 — ^(fc)] exp[— (t/r/j)^], 
the green dotted line represents the a-relaxation and satis- 
fies A{k)exp[—{t/Ta)'']- The solid green line going through 
Foo{k3,t) is the sum of both. 



essentially stuck in a cage formed by its neighbors. The 
/3-relaxation time Tp is of the order of picoseconds. On 
longer time scales, the molecule can escape from its cage 
and diffuse away from its initial position. The time Tq, is 
the relaxation time of this structural process. 

Mode-coupling theory of supercooled simple liquids 



predicts that [52] 

Foo{Kt) = Sooik, t)/Sooik,0) 

= [1 - A{k)] e-(*/^'^)' + A{k) e-^*/^")' (5) 

The factor A{k) is the Debye- Waller factor arising from 
the cage effect, which is independent of the temperature 
and follows A{k) = exp(— with a the radius of 
the cage. We are able to fit Eq. (|5| remarkably well to 
all our data, as for example in Fig. 171 

Data in Fig. [7] was collected every 10 fs for simulations 
of 1 ns. This rate of sampling results in a large amounts 
of data and is unfeasible for our runs up to 1000 ns. 
Therefore, for the 1000 ns runs we collect data at 10 ps 
intervals. At this rate of sampling it is no longer possible 
to estimate or the cage size a, but it is still possible to 
determine accurately, utilizing the fact that Soo{k,t) 
reaches a plateau near t ~ Tp. One can therefore define 

Coo{k,t) = Soo{k,t)/Soo{k,Tp), (6) 

which is S{k,t) normalized by its value at the plateau 
(Fig. [s]). A good estimate of is then the time for 
which Coo(fc,Ta) = 1/e « 0.37. 
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FIG. 8: Decay of Soo(ki,t) with time, for three different 
temperatures at P = 210 MPa {N = 343). Using Coo{k,t) 
with Tp as the approximately time when Soo{k,t) reaches a 
plateau, it is possible to obtain a good estimate of Tq. Indi- 
cated here is the Tc for 243 K, equal to ~ 60 ns. At given 
P and T we define the correlation time r as the longest time 
Ta(k) for which Coo{k,Ta) = 1/e (thin dashed black line). 

From the shorter 1 ns runs (which were mostly done 
in the HDL regime) we find that the cage radius is a = 
0.35 ± 0.09 A with a stretching exponent of 6 = 0.63 ± 
0.09. Both parameters a and b do not show a significant 
dependence on the state point within the studied range 
of temperature and pressure. 

As shown in Fig. [7] different k result in slightly differ- 
ent values for t^. We use as the correlation time r the 
largest value of which is usually found at k ^ ki, the 
first maximum in {Soo{k)) (inset Fig. [t]). 
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FIG. 9: Correlation function for four systems with different 
sizes at 210 MPa and 243 K. As the curves are all quite similar 
(which is a result we also find for the other state points), we 
conclude that the system size has a negligible effect on the 
correlation time. 



As is to be expected, the correlation time does not 
seem to depend on the box size (Fig. [9]). It does how- 
ever depend strongly on the phase, which is evident from 
Fig. [101 
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states, with correlation times in the nanosecond range. 

If we lower the temperature further, the correlation 
time slowly increases until the system becomes a glass 
rather than a liquid, and we are no longer able to fully 
equilibrate the system. As we can only run simulations 
up to 1000 ns, we consider the state points with a cor- 
relation time above 100 ns to be beyond our reach. We 
therefore designate the effective glass transition temper- 
ature Tg as the temperature for which t > 100 ns (see 
Fig.[l}. 



V. STRUCTURAL PARAMETERS 

Apart from the intermediate scattering function, there 
are other ways to quantify the structure of a liquid. In 
this section we shall examine several structural parame- 
ters, and determine which of those are the most effective 
in distinguishing between LDL, HDL, and the crystal. 
For simplicity, we approximate the center of mass of a 
water molecule with the center of its oxygen atom. 

The structural parameters are designed to distinguish 
between different phases by analyzing the geometrical 
structure. This is typically done by evaluating the spher- 
ical harmonics Y™{Lp, for a particular set of neighbor- 
ing atoms, with ip and 'd the polar angles between each 
pair of oxygen atoms in that set. In this paper we con- 
sider two different sets: we define the first coordination 
shell ni {i) to be the four nearest neighbors of molecule 
i, and define the second coordination shell n2(i) as the 
fifth to sixteenth nearest neighbors (the sixteenth nearest 
neighbors minus those in the first shell). 

Different values of i are sensitive to different symme- 
tries. The spherical harmonics with £ ^ 3, for example, 
are sensitive to a diamond structure. Those with £ = 6 
are more sensitive to the hexagonal closest packing (hep) 
structure. Since we expect the liquid and crystal struc- 
tures to be hep, diamond, or a mix of these, we focus 
primarily on £ = 3 and £ — 6. 



A. Parameters and 



FIG. 10: Arrhenius plot of the correlation time r for differ- 
ent pressures. Errors on our estimates are of the order of the 
discontinuities along the curves. At high temperatures (the 
HDL regime) the correlation time is of the order of 10-100 ps, 
which jumps several orders up as we pass the phase transi- 
tion line and enter the LDL regime. To obtain this plot, we 
dismissed the simulations that had a significant increase in r 
because of crystal growth (Sec. |VI[ ). 

At high temperatures the system is in the HDL phase, 
and has a correlation time t on the order of 10-100 ps. 
As we decrease the temperature at fixed pressure, the 
value of T has a large increase when we cross the phase 
transition line or the Widom line, depending if P is above 
or below Pc, respectively. Apparently, the LDL states 
evolve nearly four orders of magnitude slower than HDL 



All parameters defined in this section are based on 
q'fl^ii) which quantifies the local symmetry around 
molecule i. It is defined as 

E -£<m<£ (7) 

where £ and m are integers, s — 1,2 indicates the shell we 
are considering, with Ns the number of molecules within 
that shell (i.e. Ni = A for the first coordination shell, 
and = 12 for the second). is normalized accord- 
ing to / \Y™\'^ sm{'3)dipd'd = 1. We can consider qf^^ii) 
as a vector q^"-* [i) in a {M + 2)-dimensional Euclidean 
space having components Re(g|*^(i)) and lva{q^f^l^{i)) . 
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This means that we can define an inner product 

m——i 

+ Im(g(;^(^))Im(g(i(J))] (8) 

and a magnitude 



(9) 



The local parameter is one way to distinguish 

between different structures, and can be used to label 
individual molecules as LDL-like or HDL-like. We can 
convert it into a global parameter by averaging over all 
molecules, 



1e 



(10) 



In Fig. [n] we see that all global g^'*'' are sensitive to the 
difference between LDL and HDL, especially q^^^^ and 
Qq '. We conclude that the structural difference is vis- 
ible in both the first and second shell, and that LDL and 
HDL differ mostly in the amount of diamond structure of 
the first shell and the amount of hep structure in the sec- 
ond shell. This is confirmed by the histograms in Fig. [12] 
in which the largest difference between LDL and HDL is 
seen in q^^ and, next, in ql^^ ■ The latter is the param- 
eter that better discriminate with respect to the crystal 
structure. 
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200 220 240 260 280 300 
time t [ns] 

FIG. 11: Fluctuations of the density and the global struc- 
tural parameters as a function of time. The parameters are 
shown for one run using 343 molecules at 200 MPa and 248 K, 
the same as in Fig.js] Parameters q^J'\ q^\ and i/'a^'' (defined 
in the text) are as sensitive as p to the difference between 

LDL-like and HDL-like structures, while the others are more 

(2) (2) 

noisy, being Qg and much less sensitive than all the 
others. Qg"-* and iP'q \ for both s = 1 and 2, have similar 
behaviors that might be related to the temporary appearance 
of crystal-hke structures. 



B. Global parameters Qs and Qe 

An alternative approach, as used by Steinhardt et al. 
[55] . is to first average q'^/liU) over all molecules, defining 
Qe,rn = J2iLi 1i'mi'^)j then calculate the magnitude 



1/2 



E 



',m^;£,l 



(11) 



Our calculations show that the parameters Qg 



s) 



and 



Qg , with s — 1, 2, are n ot e fficient in discriminating 



between LDL and HDL (Fig. [ill), although Qq = Qg^' has 



been proposed recently as a good parameter to this goal 
[63] and consequently has been used by several authors 
[66r,68) . In particular, we observe that there is not much 

(s) 

correlation between the fluctuations of Q\ ' and those of 

the density, except for Qg^^ 

However, we confirm that Qg^'' and Qg^-* are excellent 
parameter to distinguish between the liquids (LDL and 



HDL) and the crystal, being the value of Qg*'' approxi- 
mately 10 tim es larger for the crystal than it is for the liq- 



uids (Fig. 13 1. This large increase of Qg"^ for crystal- like 



structures might be related to the few instances in Fig.[TTl 

(s) 

where an increase in Qg corresponds to a decrease of 
density (such as within interval t — 230-237 ns), con- 
sistent with the observation that the crystal-like struc- 
tures have a density comparable to the LDL structure 
and smaller than the HDL structure. 

To confirm that LDL remains a liquid in the thermody- 
namic limit, we look at how Qq changes with the system 
size. For liquids Qq scales like N^^/^ while for crys- 
tals the value Qg remains finite as TV — 00. We find 
that the probability distribution functions of QqN^^"^ 
for N = 216, 343, 512, and 729 overlap, which means 
that QqN^^'^ is independent of the system size, therefore 



Qg ~ iV"i/2 (Fig. 14). We conclude that the metastable 



LDL is not transforming into the stable crystal in the 
thermodynamic limit. This implies that the LDL and 
the crystal phase are separate by a free-energy barrier 
that is higher than kgT at the temperatures we consider 
here and that the system equilibrates to the stable (crys- 
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FIG. 12: Histograms of for £ = 3, 6 and coordination 
shells s = 1,2 at 215 MPa with iV = 343 molecules. The 
solid red (dark) curves correspond to HDL structures, and 
the solid blue (light) curves to LDL structures. The dashed 
black curve corresp ond s to the crystal structure found in run 
The parameter gg^' (a) discriminates 



VI 



C described in Sec. 
better between HDL and LDL structures, while the parameter 

(21 

Qg discriminates better between liquid-hke and crystal-hke 
structures. Parameters in (b) and (c) are much less sensitive 
to structural changes. 
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FIG. 13: Histograms of Q\ for £ — 3,6 and coordination 
shells s = 1,2 at 215 MPa with = 343 molecules. The 
symbols are as in Fig. 
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The parameter Qg*', for the first 
shell in (b) and the second in (d) , shows a clear difference be- 
tween the liquid-like structures and the crystal-Hke structure, 
but not between the two liquids. Note that scales on x-axis 
in panels (a) and (c) are one order of magnitude smaller than 
those in panels (b) and (d). As a consequence, Qa"', for the 
first shell in (a) and the second in (c), is much less sensitive 
to structural changes than Qg"'. 



tal) phase only in a time scale that is infinite with respect 
to our simulation time (1000 ns), as occur in experiments 
for metastable phases. Therefore, the LDL is a bona fide 
metastable state. Our conclusion is consistent with re- 
cent calculations by other authors 



C. Bond parameters dg and tps 



We define the bond order parameter di"'' similar to that for the crystal, most molecules have d 



the bonds within the same layer (three out of four) have 
d^;^\i,j) = — 1, while the bonds connecting atoms in dif- 



ferent layers (one out of four) have d''^\i,j 



-1/9. 



We find that the parameters d), for £ = 3, 6 and s — 1, 
2 do not distinguish well between the two different liquid- 
like structures, but that d'"^^ and dg*^ for both s = 1 
and 2 are suitable to discriminate between the crystal- 
like structure and the liquids (Fig. 15). In particular, 

0.87, and 



Ii) 



< 



defined by Ghiringhelli et al. in Rcf . 84J , where the quan- 
tity d'^\i,j) characterizes the bond between molecules i 
and j, and is designed to distinguish between a fluid and 
a diamond structure. The local parameter d^l\i,j) is 
defined as the cosine of the angle between the vectors 
ci['\i) and q:[\j): 



we therefore consider a molecule to be part of a crystal 
if at least three out of its four bonds with its nearest 



(1) 



< 



This is the same cutoff 



neighbors have d 
used by Ghiringhelli et 

The global parameter associated to d^^^ {i, j) is defined 



-0.87. 

L in [• 



as 



(12) 



1 ^ 



(13) 



with the inner product and magnitude as defined in 
Eqs. ([8]) and (|9|. 

A crystal with a perfect diamond structure has 
d'ip{i,j) = —1 for all bonds. For a graphite crystal only 



where 



is) 



(14) 
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FIG. 14: Finite size scaling of parameter Qe in the LDL 
phase (210 MPa, 243 K). The probabiUty distribution func- 
tion of QsN^^^ is independent of the system size A'', which 
means LDL scales like a liquid in the thermodynamic limit: 



is the average of d'f ^ {hj) over the first four nearest neigh- 
bors of the molecule i. We observe that each ip^/\i) has 
the same features of the corresponding with 
discriminating well be twee n the crystal-like and 

the liquids-like structures (Fig.[l6|. We observe that -^g^^ 
discrimina tes w ell between LDL-like and HDL-like struc- 



tures (Fig. 11 1, while ipg^^ for s = 1 and 2 might be able 



to emphasize the temporary appearance of crystal-like 

(s) 

structures, as noted for Ql ' . 



VI. GROWTH AND MELTING OF CRYSTAL 
NUCLEI 

In a small percentage of our simulations, the system 
was found to spontaneously crystallize. These are inter- 
esting events because spontaneous crystallization of wa- 
ter in molecular dynamics is extremely rare; only recently 
Matsumoto et al. were the first to successfully simulate 
the freezing of water on a computer |85j . Crystallization 
events in supercooled ST2 water are particularly impor- 
tant to study, as it has been proposed that LDL is un- 
stable against crystallization 

Following the discussion in Sec. |V] we define a crys- 
tal as a cluster of molecules which has three out of four 
bonds with d!ip < —0.87 and belong to the first coordi- 
nation shell of each other. In this section we shall study 
the growth and melting of these crystal nuclei, and es- 
timate the critical nucleus size needed to overcome the 
free energy barrier. The existence of this barrier allows 
us to conclude that LDL is in fact a bona fide metastable 
state with respect to the crystal. 

In Fig. [it] we show the density evolution for 11 dif- 
ferent configurations, each with 343 molecules and at 
205 MPa and 246 K. Each of these runs started at a dif- 
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Histograms of d^/^ for ^ = 3, 6 and coordination 
= 1,2 at 215 MPa with = 343 molecules. The 

in (a). 
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symbols are as in Fig 
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Apart from dg^' 



these pa- 
rameters do not distinguish well between the two different 
liquid-like structures, but dg^' and dg^' for the first shell (b) 
and the second (d) are suitable to distinguish between the 
crystal and the liquids. The parameter da'^' in (c) is remark- 
ably the same for the three structures. 



ferent initial density (between 0.85 and 0.95 g/cm'^) and 
was subsequently equilibrated to the final temperature 
and pressure using the procedure described in Sec. \n\ 
Because this state point lies close to the LLPT, we see 
phase flipping in all of them. However, the two configu- 
rations C and F display a sudden jump to a stable low 
density plateau. This is a hallmark of crystalization. We 
confirm this by calculating the size of the largest crystal 
as a function of time (Fig. 18). During most runs the 



largest crystal continuously grows and shrinks, but never 
reaches a size larger than 30 molecules. On the other 
hand, configurations C and F show a jump in crystal size 
exactly matching the jump in density. Run F ends up 
partially crystalized, while for C we find that over 90% 
of the box is crystallized in a diamond structure with a 
density of about 0.92 g/cm^ (Fig. 19). 



The correlation time increases dramatically if crystals 
appear with a size comparable to the system size, as is 
evident from Fig. [20] The correlation functions of C and 
F decay very slowly, leading to correlation times of 200- 
400 US, while the other configurations have a correlation 
time of less than 4 ns. 

For spontaneous crystallization to occur, a sufficiently 
large crystal nucleus needs to form within the liquid. Ac- 
cording to classical nucleation theory, this nucleus needs 
to reach a minimum size to prevent it from melting. We 
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FIG. 16: Histograms of i/j^"^ for £ = 3,6 and coordination 
shells s = 1,2 at 215 MPa with = 343 molecules. The 
symbols are as in Fig. 12 Each il>^ /\i) has similar features as 



the corresponding d^/\i,j) in Fig. 15 
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FIG. 17: Density vs. time near the phase transition line 
at P = 205 MPa and T = 246 K for several different con- 
figurations of A'^ = 343 molecules. This state point lies near 
the phase transition, and therefore phase flipping is seen to 
occur. Runs C and F (partially) crystallize and, at that mo- 
ment, cease to phase flip and remain stable at a low density. 



observed in many simulations that a small nucleus grovifs 
and melts, and a few runs in which the nucleus grows 
further or remains stable. Therefore, we can make an 
estimate of the critical nucleus size. 

The two largest crystals that formed and subsequently 
melted, both reached a size of about 50-60 molecules 
(Fig. 21 1 and 21 d). The smallest crystal that formed 
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FIG. 18: Evolution of crystal size with time for the same 
configurations as in Fig. |17[ The j/-axis goes from to 34, 
except for configurations C and F which go up to 343. The 
system spontaneously crystallizes in both C and F, while the 
largest crystals in the remaining configurations never reach a 
size larger than 30 molecules. 




FIG. 19: A snapshot (at t = 1000 ns) of the cubic diamond 
crystal produced by run C of Figs. [l7] and |18| Shown here 
are all A'' — 343 molecules, with a small part still in the liquid 
state (bottom-left corner), and a crystal defect in the center. 
Note that the defect only affects the position of the hydrogen 
atoms, and not that of the oxygen. 



and remained stable, had a size of about 50-80 molecules 
(Fig. 21:). We therefore conclude that the critical nu- 
cleus size is approximately 70 ± 10 molecules. A similar 
value of ~ 85 molecules was found by Reinhardt and Doye 



for ice nucleation in the monatomic water model [57] . 

For a more accurate estimate it is necessary to run 
longer simulations, as the crystal nuclei can survive for 
hundreds of nanoseconds (e.g.. Fig. 21 i in which a small 
crystal lasts for 700 ns). 
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FIG. 20: The correlation time increases dramatically if crys- 
tals of a size com para ble t o th e system size appear (i.e. runs 
C and F of Figs. 17 and 1181). The correlation time of two 



other runs (H and Jj are slightly larger than average because 
these runs spend more time in the LDL phase (see Fig. |17|l . 
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FIG. 21: Growth and melting of crystal nuclei, (a) The 
largest nucleus that melted reached a size of 62 molecules 
during a simulation of 512 molecules at 210 MPa and 244 K. 
(b) The second-largest nucleus was 55 molecules during a sim- 
ulation of 343 molecules at 210 MPa and 243 K. (c) A few 
runs lead to irreversible crystallization (A*' = 216 at 195 MPa 
and 245 K). (d) Some crystal nuclei survive for hundreds of 
nanoseconds (A'^ = 343 at 195 MPa and 246 K) before disap- 
pearing. 



VII. LOCATING THE CRITICAL POINT 



In Sec. |III| we used the intermediate scattering func- 
tion Sooik) to estimate the position of the liquid- liquid 
critical point, and found it to lie near 200-210 MPa and 
244-247 K. It is commonly believed that the LLCP falls 
in the same universality class as the three-dimensional 
Ising model [57] . At the critical point the order param- 
eter distribution function (OPDF) of a system has the 



same bimodal shape as all other systems that belong to 
the same universality class. Therefore we can locate the 
LLCP accurately by fitting our data to the OPDF of the 
3D Ising model (Fig. [22]). 
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FIG. 22: Fitting the order parameter distribution function 
(obtained from the simulation data) to that of the 3D Ising 
model at the critical point (black curve, from Ref. 88 ). The 
shape of the order parameter distribution function depends 
on temperature and pressure. Here (for N = 343) we find 
an excellent fit for Pc = 206 MPa and Tc = 246 K, and 
we therefore confirm that the LLCP indeed belongs to the 
same universality class as the 3D Ising model. Based on our 
fit, we locate the LLCP to be at Pc = 206 ± 3 MPa and 
Tc = 246 ± 1 K. 

In the 3D Ising model the order parameter M is simply 
the spontaneous magnetization, but for liquids the order 
parameter turns out to be a linear combination of two in- 
dependent quantities such as the density and the poten- 
tial energy [89, 90 . We therefore define M = p+sE, with 
s a constant known as the field mixing parameter. The 
value of s depends only on the model and should therefore 
be independent of the number of molecules. Our fits in- 
deed confirm this; we find s = 0.0362 (g/cm^)/(kJ/mol) 
for all values of N. 

Only the shape of the OPDF is dictated by the theory, 
which means we are free to move and stretch our OPDF 
to acquire an accurate fit. So, instead of fitting the order 
parameter M, we actually fit a; = A{M — Mc) to the 3D 
Ising model. The critical order parameter Mc is chosen 
such that the mean value of x is zero, and the amplitude 
A has been chosen such that the variance equals unity. 

To calculate the OPDF for a particular pressure and 
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temperature, we create a two-dimensional histogram of 
the density and energy. An example of such a 2D his- 
togram is shown in Fig. [23] for P = 200 MPa, T = 
247.5 K, and iV = 343. Near the critical point the his- 
togram displays two peaks, one for LDL and one for HDL. 
If we integrate the 2D histogram along the direction cor- 
responding to the value of s, we obtain the histogram for 
M. 
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FIG. 23: 2D histogram of the density and the total energy 
for a system at 247.5 K and 200 MPa (on the Widom line), 
obtained via histogram reweighting. The histogram of the 
energy (curve a) seems to indicate that the system is mostly 
in the LDL state, while the histogram of the density (curve 
b) indicates the HDL state is more predominant. For liquids 
the order parameter M = p -\- sE is actually a linear com- 
bination of the density p and the energy E (curve c), with 
s = 0.0362 (g/cm3)/(kJ/mol). 

To fit our OPDF to that of the 3D Ising model, we need 
to calculate our OPDF at different pressures and temper- 
atures, until we find the {Pc,Tc) that gives us the best 
fit. The state point {Pc, Tc) is then our best estimate of 
the location of the LLCP. We only have simulation data 
for a finite number of state points, therefore some kind 
of interpolation is necessary. The method of choice here 
is histogram reweighting ^T\; we use the algorithm as 
described by Panagiotopoulos in Ref. [92] , 

The results of fitting our data to the 3D Ising model 
are shown in Fig. [24j Tables |T] and [IT] indicate which 
data was used by the histogram reweighting method to 
obtain these fits. For N — 343, 512, and 729, we are 
able to fit our data very accurately to the OPDF of the 
3D Ising model, and find the critical point to be located 
at Tc = 246 ± 1 K, Pc = 206 ± 3 MPa for N = 343, 
and at Tc = 246 ± 1 K, Pc = 208 ± 3 MPa for N ^ 
529 and 729. Theory predicts that the location of the 
critical point depends on N, and these findings agree 
with that prediction. In particular, the 3D Ising model 
predicts that the amplitude A should scale with box size 
L as A ~ L^/"" oc N^^^" with (3/i^ = 0.52 [Ml ESI, in 
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FIG. 24: Best OPDF fits for iV = 343, 512, 

and 729 molecules. In all cases we obtain s = 
0.0362 (g/cm'^)/(kJ/mol) for the field mixing parameter. We 
find the critical point to be located at Tc ~ 246 ± 1 K, 
Pc = 206 ± 3 MPa for TV = 343, and at Tc = 246 ± 1 K, 
Pc = 208 ± 3 MPa for N = 529 and 729. 
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TABLE I: Number of simulations used for the histogram 
reweighting in order to obtain the order parameter distribu- 
tion function for TV = 343 molecules. 



also indicates that TV — 216 cannot provide an accurate 
estimate of the location of the LLCP. 

To establish that the LLPT does not vanish in the 
thermodynamic limit TV — >■ (X), we consider the finite size 
scaling of the Challa-Landau-Binder parameter [94U99) . 
Near the critical point the density distribution function 
T>{p) has a bimodal shape that can be approximated by 
the superposition of two Gaussians (e.g.. Fig. 23). The 



Challa-Landau-Binder parameter 11 is a measure of the 
bimodality of T>{p) and is defined as 



n = 1 - 



3{pr 



(15) 



agreement with the slope of A(N) in Fig. 25 This figure 



When there is only one phase, T>{p) is unimodal and 11 = 
2/3. But in a two-phase region, with two phases that 
have different densities, the shape of T>{p) is bimodal 
(Fig.[23| and 11 < 2/3. For a finite system T>{p) is always 
bimodal at both the Widom line and the LLPT, but in 
the thermodynamic limit there exists only one phase at 
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TABLE II: Simulations used for the histogram reweighting 
in order to obtain the order parameter distribution function 
for N = 512 and N = 729 molecules. 
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FIG. 25: Log-log plot of the amplitude A vs. system size N. 
From the slope of this line we determine that A ~ A'''^'^^ oc 
1/° **, in fare agreement with the value of 0.52 predicted by the 
3D Ising model :89 . For the smaller size A'' = 216 we observe 
large finite-size deviation with respect to the thermodynamic 
limit behavior. 



the Widom line, while there remain two at the phase 
transition hne. Therefore, 11 — >■ 2/3 at the Widom hne, 
while n < 2/3 at the LLPT even in the limit N oo. 
Hence, the finite-size scaling of 11 allows us to distinguish 
whether an isobar crosses the LLPT or the Widom hne, 
and is yet another method of estimating the location of 
the critical point. 

We study 11 versus temperature T and system size N 
for different pressures, finding minima IXmin at specific 
temperature for each pressure (Fig, 
dependence of nnii„(P) reveals if P 



261 



< 



mm 

The finite-size 
Pc or P > Pc 



(Fig. 271 



For P < Pc the mimimum Umin approaches 2/3 hn- 
early with 1/A'', while for P < Pc it approaches the limit 




235 240 
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FIG. 26: The Challa-Landau-Binder parameter 11 as a func- 
tion of temperature and system size A^, for four different pres- 
sures. For finite system sizes 11 shows a minimum at the 
LLPT and the Widom line, while 11 2/3 (thin dashed hne) 
at temperatures where ©(p) is given by a single Gaussian. 
The finite-size scaling of the minimum of II, indicates that 
the critical point exists in the thermodynamic limit (Fig. 27 1 . 



2 ^2 



1 ipii-pj) 



(16) 



This limiting value is also approached linearly with 1/iV. 
Here pn = Pn{P) and pl = Pl{P) are the densities of 
the two phases LDL and HDL [9W. Above the critical 
pressure the limiting value of Umin decreases as P in- 
creases, i.e. the two peaks of the bimodal ^{p) move 
further apart. This happens because pu ~ Pl increases 
at coexistence as {P — Pc)^ where /3 w 0.3 is the critical 
exponent of the 3D Ising universality class [1001 llOlj . 

From this analysis (Fig. 27 1 we conclude that our re- 
sults agree with theory and that the critical pressure 
Pc « 190-210 MPa, in agreement with the estimate of 
Sec. VII Furthermore, as H remains less than 2/3 for 
P > Pc even in the limit N oo, we conclude that the 
LLPT does not vanish in the thermodynamic limit. 



VIII. CONCLUSIONS 

We performed molecular dynamic simulations in the 
NPT ensemble for ST2-RF water in the supercooled re- 
gion of the phase diagram for different system sizes with 
simulation times of up to 1000 ns. Using several differ- 
ent techniques we confirmed the existence of two liquid 
phases, LDL and HDL, separated by a liquid-liquid phase 
transition line. Near the LLPT line the system continu- 
ously flips between the two phases. Because of finite size 
effects this phenomenon also occurs near the Widom line, 
but by fitting the order parameter distribution function 
to that of the 3D Ising model, we were able to accurately 
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0.654, 
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FIG. 27: Minima of the Challa-Landau-Binder parameter IT 
as a function of system size for different pressures. The 
minimum Ilmin occurs at the pressures and temperatures of 
the LLPT and the Widom line, and is always less than 2/3 
for a finite system because of the bimodality of the density 
histogram. As — >■ oo the bimodality disappears in the one- 
phase region but remains at the LLPT, and therefore flmin — >■ 
2/3 at the Widom line while Ilmin < 2/3 on the LLPT, even in 
the thermodynamic limit. We conclude that the critical point 
survives in the thermodynamic limit, and that it is located 
between P — 200 and 210 MPa (in agreement with previous 
results of this paper). 



determine the location of the liquid-liquid critical point 
(at Tc = 246 ± 1 K, Pc = 208 ± 3 MPa). Finite size 
scaling of the Challa-Landau-Binder parameter indicates 
that the critical point docs not disappear in the thermo- 
dynamic limit. 

Both phases have been confirmed to be bona fide 
metastable liquids that differ substantially in structural 
as well as dynamical properties. It is found that the LDL 
phase is a more "structured" liquid, and that it has a cor- 
relation time of almost four orders of magnitude larger 



than that of HDL, with LDL correlation time of the order 
of 100-1000 ns. We show that Qe structural parameter is 
not able to discriminate between HDL and LDL, but can 
discriminate well between liquids and crystal. Finite size 
scaling of the Qq parameter confirms that LDL scales as 
a liquid and not as a crystal. 

The different structures of LDL and HDL are better 
discriminated by structural parameters like q^^^^ and ql^^ ■ 
These parameters show that LDL and HDL differ mostly 
in the amount of diamond structure of the first shell and 
the amount of hep structure in the second shell. 

For small box sizes {N = 343) there were a few simu- 
lation runs that resulted in spontaneous crystallization, 
always within the LDL region of the phase diagram. Fur- 
ther analysis revealed that during all simulations small 
crystals grow and melt within the liquid, a clear indica- 
tion that LDL is metastable with respect to the crystal. 
From the few crystalization events that occurred, we were 
able to conclude that the critical nucleus size is approxi- 
mately 70 ± 10 molecules. 
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